Open shells in reduced-density-matrix-functional theory 
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Reduced-density-matrix-functional theory is applied to open-shell systems. We introduce a spin- 
restricted formulation by appropriately expressing approximate correlation-energy functionals in 
terms of spin-dependent occupation numbers and spin- independent natural orbitals. We demon- 
strate that the additional constraint of total-spin conservation is indispensable for the proper treat- 
ment of open-shell systems. The formalism is applied to the first-row open-shell atoms. The obtained 
ground-state energies are in very good agreement with the exact values as well as other state of the 
art quantum chemistry calculations. 
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The proper description of electronic correlation is the 
central problem in theoretical material science. Den- 
sity functional theory (DFT) deals with this problem 
by considering the electronic density as the fundamen- 
tal variable. DFT is built upon the Hohenberg-Kohn 
theorem 1] which ensures that the expectation value of 
every observable, including the total energy, is a unique 
functional of the electronic ground-state density. The 
many-electron problem is then mapped onto an auxil- 
iary, non-interacting system, the so-called Kohn-Sham 
system^ whose density reproduces the interacting den- 
sity. In practice, the so-called exchange and correlation 
part of the Kohn-Sham potential needs to be approx- 
imated. DFT is tremendously successful in describing 
properties of real materials and that success is proved 
by its general acceptance as one of the major tools in 
exploring theoretically the world of real materials. Nev- 
ertheless, there are certain materials, termed collectively 
as highly correlated materials, for which the results of 
DFT, at least with standard functionals for exchange and 
correlation, deviate significantly from the experimental 
values. 

Reduced-density-matrix-functional theory (RDMFT) 
is an alternative way to deal with the many-electron prob- 
lem. It is based on Gilbert's theoremyf, which guaran- 
tees that the expectation value of any observable of the 
system in the ground-state is a unique functional of the 
one-body reduced-density-matrix (1-RDM) 

7 (r, r') = Njd 3 r 2 - ■ -d 3 r N ^*(r', r 2 - • -r N )V(r, r 2 - • -r N ) 

(1) 

where VP is the many-body ground-state wavefunction. 
The advantage of this approach, compared to DFT, is 
that the exact many-body kinetic energy is easily ex- 
pressed in terms of 7. Although the properties of 1-RDM 
functionals had been the subject of theoretical studies for 
a long time0, 0, 0, 0, there were few practical applica- 
tions of RDMFT until recently. 

We emphasize that it is not possible to con- 
struct a Kohn-Sham-like independent electron scheme in 
RDMFT: Due to its non-idempotency, the 1-RDM of 



an interacting system cannot be reproduced by a non- 
interacting system because the latter always has an idem- 
potent 1-RDM. This is reflected in the eigenvalues of 
the 1-RDM which are equal to zero or one for a non- 
interacting system, while in the interacting case some of 
them arc fractional. 

The first implicit functional of the 1-RDM was intro- 
duced in 1984 by Miiller 6]. This functional was later 
derived from modeling the correlation hole, by Buijse 
and Baerendsjj]]- Goedecker and Umrigar|3| introduced a 
self-interaction correction to this functional. Other func- 
tionals have been proposed more recently^ H [H E 
1 1.1 Il4j . Implicit 1-RDM functionals depend explicitly on 
the so-called natural orbitals <p a and the corresponding 
occupation numbers n a which are defined as the eigen- 
functions and the eigenvalues of the 1-RDM 



dr' 3 7(r, r') ip a (r') = n a ip a (r) 



(2) 



The total energy can be minimized with respect to the 
natural orbitals and the occupation numbers, instead 
of the 1-RDM itself, by considering the appropriate 
subsidiary conditions, namely the conservation of the 
total number of electrons, the JV-representability con- 
straint and the orbital orthonormality. Goedecker and 
UmrigarQ calculated the correlation energy of small 
atomic systems including the open-shell Lithium and 
Carbon atom. Despite their simplicity, 1-RDM function- 
als are able to provide a good approximation of the corre- 
lation energy of small systems H0HHE1E1E1E1- 
Generally, the correlation energy of finite systems is im- 
proved by removing the self- interaction terms fl5t fl6| . 
However the correct dissociation limit of the Ho molecule 
is found only if self-interaction is retained^. Il5l . 

In the present work, we demonstrate how natural or- 
bital functionals in RDMFT can be applied to open-shell 
systems. For that purpose, we employ the Miiller func- 
tional @ with the self-interaction correction proposed 
by Goedecker and Umrigar||. Our formulation is spin- 
restricted, i.e. we have spin-dependent occupation num- 
bers but spin-independent natural orbitals. The advan- 
tage of this approach is that particular spin configu- 
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rations are prescribed, in the same manner as in the 
restricted open-shell Hartree-Fock (ROHF) method^}- 
We discuss the necessity of a spin-dependent constraint 
for the conservation of both the total number of electrons 
as well as the spin. 

We start with the 1-RDM functional described in 
Refs. 0, 0, 0|> where we explicitly include the full spin 
dependence in the occupation numbers and the natural 
orbitals. We refer to this as the spin-unrestricted func- 
tional that reads 



E 



\ «'(i-<u«W)J 7' 

ab,a<j' 

- \Y.VK^^~^ b )K b , (3) 



ab.a 



with being the spin-dependent occupation numbers 
and h^ a , J^b an d K^ b the one- and two-electron integrals 
for the natural orbitals (p aa 



V 2 

-— + V ext {v) 



</W( r ) , 



J ab 



* d 3 r jy ^LW y ag (r) yw(r') <p* ba ,(r') ^ 



In the above expression, V ex t (r) is the external potential, 
i.e. the ionic potential for atomic and molecular systems. 
The Kronecker Ss are inserted in Eq. in order to ex- 
clude explicitly the self-interaction terms 0- 

We now concentrate on the spin-restricted case 
by assuming spin-independent orbitals but still spin- 
dependent occupation numbers. Then the expression @ 
reduces to 

E = ( n i + U a) h ™+ 

a 

\ + n i n b) f 1 - S ab) + (n ] a n l b + n l a nl^j J ab 



-^2 (\j n * n i + \j 7l « ni b) (1-Sab) 



ab 



K ab . (5) 



In this equation, the integrals h aa , J a b and K ab are analo- 
gous to those in Eq. Q but for spin-independent orbitals 
ip a . This is the functional we employ in this work. It can 
be viewed as a generalization of the ROHF functional in 
the sense that it reduces to the ROHF functional if all the 
occupation numbers are either zero or one. Although the 
formulation in this paper refers to the specific functional 
given in Eq. © , the generalization to different function- 
als of similar kind[2j is easily achieved by replacing the 
square roots in the exchange-like term in Eq. <0 by the 
appropriate functions f{n\,n b ), f(n a ,n b ). 



Driven by physical requirements, Cioslowski et al .18] 
have derived a list of criteria that 1-RDM functionals 
should fulfill. In addition, they examined whether func- 
tionals of the Miiller-type satisfy these criteria. Our 
open-shell version JSJ of the self-interaction-corrected 
Miiller functional satisfies the same criteria as the original 
Miiller functional^ with the self-interaction correction of 
Gocdccker and Umrigar 8]. 

We now discuss the subsidiary conditions that have to 
be enforced in the minimization procedure of function- 
als like the one given by Eq. (JSJ). Since the occupation 
numbers are spin-dependent we face a dilemma concern- 
ing the conservation of the total spin. More specifically 
we have the following two options: (i) to use one con- 
straint for the conservation of the total number of elec- 
trons. This introduces a single Lagrange multiplier, the 
chemical potential /z. In that way, we unfortunately al- 
low for charge transfer from one spin channel to the other, 
usually from the majority to the minority spin and the 
total spin is not preserved, or (ii) to use an extra con- 
straint for the conservation of the total spin, making the 
minimization more restrictive. In practice, we use two 
different constraints for the spin-up and the spin-down 
electrons. In that way, we introduce two Lagrange mul- 
tipliers, or, in other words, a spin-dependent chemical 
potential. It is one of the goals of the present work to 
assess these two different ways of minimizing the energy 
functional. Of course, the above dilemma applies only 
to open-shell systems. It is expected that the first op- 
tion, being less restrictive, leads to a lower total energy. 
However, since RDMFT (like DFT) is not variational, a 
lower energy is not necessarily better. Hence, it is not a 
priori clear whether enforcing the spin conservation con- 
straint in addition to the particle number conservation 
will improve or worsen the energy. 

Additionally, as in the case of closed-shell systems, we 
have to include two further subsidiary conditions. The 
first is the orbital orthonormality. Unfortunately, all 1- 
RDM functionals are not invariant under unitary trans- 
formation of the natural orbitals. This leads to a complex 
minimization problem which consists of orbital depen- 
dent Fock-like operators and off-diagonal Lagrange mul- 
tipliers. The equations we have to solve in order to find 
the orbital solution for fixed but fractional occupation 
numbers have the form 

F< >(r)p«(r) = ]Te Qb ^(r), 



^ba — tab ; 



(6) 



where is analogous to the Fock matrix in Hartree- 
Fock theory but in this case orbital- index dependent, 
and tab are the non-diagonal Lagrange multipliers, which 
should be Hermitian according to the second equation. 
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The orbital-index dependent operator is 
F^(r) = (nl+ni)h(r)+ 

[( U a n l + nl a n i) (1 - S ab) + (™M + flj;nj)] Jb(r) 



XI v n « n fc T + v nlaU b ^ ~ 5ah ) 



where 



ab 

h(r) 
Mr) 
K. b {r)(fa(r) 



Mr), (7) 



V 2 



rfr' 



r - r 



(8) 



r - r' 



A similar problem with an orbital-index dependent oper- 
ator exists in the self-interaction-correction met hod (ill 
Hj of DFT and to a lesser extent in the ROHF 
methodfl?! Ell ]. Finally, the second of the subsidiary 
conditions, known as the N-representability constraint, 
restricts the occupation numbers to lie between zero and 
one, < n a a < 1 and guarantees that the 1-RDM corre- 
sponds to an iV-body wave function. 

Our implementation uses the GAMESS quantum 
chemistry program 23J for the calculation of the one- 
body and two-body integrals as well as for providing the 
Hartree-Fock solution which we choose as the starting 
point of our minimization procedure. The minimizations 
with respect to both the occupation numbers and the 
natural orbitals are performed with a conjugate gradi- 
ent procedure. For the optimization with respect to the 
orbitals, we ado pte d a procedure similar to the one de- 
scribed in Refs. [2fJ,|23. We calculated the total energies 
of the first-row open-shell atoms, i.e. Li, B, C, N, O, and 
F. For all these elements we used the cc-PVQZ Gaussian- 
basis set [24). 

In Figure ^ we show the convergence of the total en- 
ergy as a function of the number of natural orbitals in- 
cluded in the calculation. As it can be seen, 30-40 natural 
orbitals are typically required for full convergence. In Ta- 
ble J] we list the total energies for the open-shell systems 
as well as the He and Be closed-shell atoms for complete- 
ness. Comparing the total energies it becomes apparent 
that forcing or relaxing the constraint of spin conserva- 
tion is extremely important for open-shell systems. The 
energy differences between the two minimization proce- 
dures are of the order of 10 mH with the spin-conserving 
constraint giving results much closer to the exact ones. 
Relaxing the constraint results in a charge transfer from 
one spin to the other which is of the order of 0.05-0.1 elec- 
trons. This charge transfer increases with the number of 
natural orbitals included in the minimization procedure 
which is the reason for the increase in the energy dif- 
ference between the two procedures with the number of 
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FIG. 1: The convergence of the total energy with the number 
of natural orbitals included in the minimization procedure for 
Lithium and Carbon atoms. The values of the ROHF method 
as well as the CI and exact values are shown as horizontal 
lines. 



naturals orbitals seen in Fig.^ Interestingly, for the last 
two open-shell elements, i.e. O and F, where the 2p shell 
is more than half filled, the energy difference between the 
two minimization procedures is smallest. 

It is clear from Fig. □ and Tabic [Q that RDMFT offers 
a very good approximation for the correlation energy for 
all, closed and open-shell, atomic systems we studied. 
Despite its simplicity, the functional we used produces 
results in very good agreement with state of the art, but 
computationally much more expensive, methods like CI. 
However, the systematic trend of this particular func- 
tional is to overestimate slightly the correlation energy 
for all the open-shell systems studied. To give a fair credit 
to RDMFT one has to take into account that this is one 
of the very first 1-RDM functionals that has been ex- 
ploited in a minimization procedure, in contrast to DFT 
functionals that have been heavily used and tuned for 
decades. 

In Table |n] we quote the different values of the chem- 
ical potential for the majority and the minority spin as 
well as the common value of it in the case of relaxing 
the constraint of spin conservation. Clearly, the values 
are significantly different. Interestingly, for B and C with 
the p-shell less than half filled, but also for N with half 
filled p-shell, the majority fx' is larger than the minority 
fi^ . This is a consequence of the fact that ionization from 
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Atom 


RDMFT 


QCI 


Exact 




T I 

fj,' = n* 


fi l / // 






He 


2.8978 


2.9049 


2.9037 


Li 


7.4940 


7.4827 


7.4743 


7.4781 


Be 


14.6686 


14.6657 


14.6674 


B 


24.6746 


24.6616 


24.6515 


24.6539 


C 


37.8675 


37.8506 


37.8421 


37.8450 


N 


54.6096 


54.5965 


54.5854 


54.5893 





75.0668 


75.0671 


75.0613 


75.0670 


F 


99.6951 


99.6952 


99.727 


99.734 


A 


0.020 


0.010 


0.004 


0.000 



good agreement with state-of-the-art quantum-chemistry 
calculations. The presented formalism can be easily ex- 
tended to any 1-RDM functional and therefore opens up 
the wide range of open-shell systems to the testing of new 
1-RDM functional. 

We thank A. Zacarias for helpful suggestions. This 
work was supported in part by the Deutsche Forschungs- 
gemeischaft within the program SPP 1145, by the EX- 
CITING Research and Training Network and by the 
NANOQUANTA Network of Excellence. 



TABLE I: Absolute total energies for the first row atoms 
(in a.u.) with and without enforcing the spin conservation 
constraint. The and exact 26] values are also listed. 

A is the mean absolute deviation from the exact values (for 
open-shell systems only). 



TABLE II: The values of the spin-dependent chemical poten- 
tial /x T (the majority spin) and ^ in the case of enforcing the 
spin conservation constraint, and the common value of jj, and 
the charge transfer Aq in the case of relaxing that constraint. 



the majority is energetically favorable. For the remain- 
ing two elements, i.e. O and F, with more than half filled 
p-shell this trend is opposite. 

In summary, we have presented a systematic appli- 
cation of RDMFT to open-shell systems. We adopted 
a spin-restricted open-shell treatment and extended 1- 
RDM functionals to include spin-dependent occupation 
numbers. This formalism has the advantage of allowing 
the prescription of a specific spin state. In particular, we 
introduced a spin-dependent chemical potential in order 
to enforce conservation of the total spin in the minimiza- 
tion procedure, in addition to the conservation of the 
total number of electrons. We performed calculations 
for open-shell atoms using the appropriate extension of 
a standard 1-RDM functional [(Jig, both with and with- 
out the spin conservation constraint. We conclude that 
the inclusion of this constraint is essential for the proper 
treatment of open-shell systems. Our results for the to- 
tal energies of the first row open-shell atoms are in very 
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Ag [e] 


[6] 
[7] 


He 


0.928 


0.928 


0.928 




Li 


0.191 


2.541 


0.186 


0.012 




Be 


0.292 


0.292 


0.292 




[8] 


B 


0.234 


0.421 


0.193 


0.089 


[9] 
[10] 


C 


0.312 


0.550 


0.286 


0.074 


N 


0.445 


0.710 


0.411 


0.058 


O 


0.474 


0.404 


0.418 


0.075 


[11] 


F 


0.570 


0.527 


0.510 


0.060 


[12] 
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